# -*- coding: UTF-8 -*-
from abaqus import *
from abaqusConstants import *
import regionToolset

def TTmodelling(T1_pair,T2_pair,T3_pair,length,mdbname):
    '''创建三个圆管的实体模型,并对称地切分
    
    参数：
    T1_pair:以元组(直径、壁厚)格式输入，
            例如(20,1)表示生成直径20，壁厚1的圆管。
    T2_pair和T3_pair同理
    length:以数字输入，表示各管的长度
    mdbname:以字符串输入，是模型数据库的名字
    
    '''

    #======创建Model=========
    currModel = mdb.Model(name= mdbname)
    #=================创建Part===============
    T1Sketch = currModel.ConstrainedSketch(name='__profile__', sheetSize=200.0)
    T1Sketch.CircleByCenterPerimeter(center=(0,0), point1=(0,T1_pair[0]/2.0))
    T1Sketch.CircleByCenterPerimeter(center=(0,0), point1=(0,T1_pair[0]/2.0-T1_pair[1]))
    T1Part = currModel.Part(name='T1' , dimensionality = THREE_D , type = DEFORMABLE_BODY)
    T1Part.BaseSolidExtrude(sketch=T1Sketch, depth=length)
    del currModel.sketches['__profile__']
    T1Part.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=0.0)
    T1Part.DatumPlaneByPrincipalPlane(principalPlane=XZPLANE, offset=0.0)
    T1Part.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=length/2)
    # 基准平面是字典
    datums = T1Part.datums.values()
    for datum in datums:
        T1Part.PartitionCellByDatumPlane(datumPlane=datum, cells=T1Part.cells)

    T2Sketch = currModel.ConstrainedSketch(name='__profile__', sheetSize=200.0)
    T2Sketch.CircleByCenterPerimeter(center=(0,0), point1=(0,T2_pair[0]/2.0))
    T2Sketch.CircleByCenterPerimeter(center=(0,0), point1=(0,T2_pair[0]/2.0-T2_pair[1]))
    T2Part = currModel.Part(name='T2', dimensionality=THREE_D, type=DEFORMABLE_BODY)
    T2Part.BaseSolidExtrude(sketch=T2Sketch, depth=length)
    del currModel.sketches['__profile__']
    T2Part.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=0.0)
    T2Part.DatumPlaneByPrincipalPlane(principalPlane=XZPLANE, offset=0.0)
    T2Part.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=length/2)
    datums = T2Part.datums.values()
    for datum in datums:
        T2Part.PartitionCellByDatumPlane(datumPlane=datum, cells=T2Part.cells)

    T3Sketch = currModel.ConstrainedSketch(name='__profile__', sheetSize=200.0)
    T3Sketch.CircleByCenterPerimeter(center=(0,0), point1=(0,T3_pair[0]/2.0))
    T3Sketch.CircleByCenterPerimeter(center=(0,0), point1=(0,T3_pair[0]/2.0-T3_pair[1]))
    T3Part = currModel.Part(name='T3', dimensionality=THREE_D, type=DEFORMABLE_BODY)
    T3Part.BaseSolidExtrude(sketch=T3Sketch, depth=length)
    del currModel.sketches['__profile__']
    T3Part.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=0.0)
    T3Part.DatumPlaneByPrincipalPlane(principalPlane=XZPLANE, offset=0.0)
    T3Part.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=length/2)
    datums = T3Part.datums.values()
    for datum in datums:
        T3Part.PartitionCellByDatumPlane(datumPlane=datum, cells=T3Part.cells)

def analyticalPlate(rectData,mdbname):
    '''生成一个矩形的刚性压板，并且在中心设置参考点
    
    参数：
    mdbname:字符串，模型数据库的名字
    rectData:元组，(长，宽)
    
    '''
    
    currModel = mdb.models[mdbname]
    plateSketch = currModel.ConstrainedSketch(name='__profile__', sheetSize=200.0)
    halfLen = rectData[0]/2.0
    plateSketch.Line(point1=(-1.0*halfLen, 0.0), point2=(1.0*halfLen, 0.0))
    platePart = currModel.Part(name='rigidPlate', dimensionality=THREE_D, 
                                                type=ANALYTIC_RIGID_SURFACE)
    platePart.AnalyticRigidSurfExtrude(sketch=plateSketch, depth=rectData[1]*1.0)
    del currModel.sketches['__profile__']
    platePart.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=0.0)
    platePart.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=0.0)
    # 基准平面是字典
    datums = platePart.datums.values()
    for datum in datums:
        platePart.PartitionFaceByDatumPlane(datumPlane=datum, faces=platePart.faces)
    platePart.ReferencePoint(point = (0.0,0.0,0.0))

def TCNTAssembly(T1_pair,T2_pair,T3_pair,tube_len,rectData,mdbname):
    '''将三个管的装配为TCNT形式, 
    并且在分析刚体上建立两个Surface:Surf_topPlate, Surf-bottomPlate
    和分别包含其上的参考点的两个Set:RF-topPlate,RF-bottomPlate
    
    接触区域绑定处理

    调用了TTmodelling()和analyticalPlate()
    
    参数：
    T1_pair:以元组(直径、壁厚)格式输入，
            例如(20,1)表示生成直径20，壁厚1的圆管。
    T2_pair和T3_pair同理
    rectData:元组，(长，宽)    
    mdbname:字符串，模型数据库的名字
    '''
    TTmodelling(T1_pair,T2_pair,T3_pair,tube_len,mdbname)
    analyticalPlate(rectData,mdbname)
    
    currModel = mdb.models[mdbname]
    currAssembly = currModel.rootAssembly
    
    #=========导入部件==================
    currAssembly.DatumCsysByDefault(CARTESIAN) # 附带默认坐标系
    T1Instance = currAssembly.Instance(name='T1-1', part=currModel.parts['T1'],dependent=OFF)
    T2Instance = currAssembly.Instance(name='T2-1', part=currModel.parts['T2'],dependent=OFF)
    T3Instance = currAssembly.Instance(name='T3-1', part=currModel.parts['T3'],dependent=OFF)
    topPlateInstance = currAssembly.Instance(name='rigidPlate-1', part=currModel.parts['rigidPlate'],dependent=ON)
    bottomPlateInstance = currAssembly.Instance(name='rigidPlate-2', part=currModel.parts['rigidPlate'],dependent=ON)
    
    #===============创建参考点SET和指定二维板的表面=====================
    r1 = currAssembly.instances['rigidPlate-1'].referencePoints
    refPoints1=(r1.values(), ) # 参考点是字典
    currAssembly.Set(referencePoints=refPoints1, name='RF-topPlate')
    r2 = currAssembly.instances['rigidPlate-2'].referencePoints
    refPoints2=(r2.values(), ) # 参考点是字典
    currAssembly.Set(referencePoints=refPoints2, name='RF-bottomPlate')
    
    s1 = currAssembly.instances['rigidPlate-1'].faces
    halflen0 = rectData[0]/4.0
    halflen1 = rectData[1]/4.0
    topPlate_surf = s1.findAt(((halflen0,-1e-7,halflen1),),((-1.0*halflen0,-1e-7,halflen1),),((halflen0,-1e-7,-1.0*halflen1),),((-1.0*halflen0,-1e-7,-1.0*halflen1),),)
    currAssembly.Surface(side2Faces=topPlate_surf, name='Surf-topPlate')
    
    s2 = currAssembly.instances['rigidPlate-2'].faces
    bottomPlate_surf = s2.findAt(((halflen0,1e-7,halflen1),),((-1.0*halflen0,1e-7,halflen1),),((halflen0,1e-7,-1.0*halflen1),),((-1.0*halflen0,1e-7,-1.0*halflen1),),)
    currAssembly.Surface(side1Faces=bottomPlate_surf, name='Surf-bottomPlate')
    
    #============找到施加Tie约束的结点域=========================
    T1_raduis = T1_pair[0]/2.0
    T2_raduis = T2_pair[0]/2.0
    T3_raduis = T3_pair[0]/2.0
    
    edge_T1_inner_top = T1Instance.edges.findAt(((0,T1_raduis-T1_pair[1],tube_len*1.0/4.0),),((0,T1_raduis-T1_pair[1],tube_len*3.0/4.0),))
    edge_T1_inner_bottom = T1Instance.edges.findAt(((0,-abs(T1_raduis-T1_pair[1]),tube_len*1.0/4.0),),((0,-abs(T1_raduis-T1_pair[1]),tube_len*3.0/4.0),))
    
    edge_T2_outter_top = T2Instance.edges.findAt(((0,T2_raduis,tube_len*1.0/4.0),),((0,T2_raduis,tube_len*3.0/4.0),))
    edge_T2_outter_bottom = T2Instance.edges.findAt(((0,-abs(T2_raduis),tube_len*1.0/4.0),),((0,-abs(T2_raduis),tube_len*3.0/4.0),))
    
    edge_T3_outter_top = T3Instance.edges.findAt(((0,T3_raduis,tube_len*1.0/4.0),),((0,T3_raduis,tube_len*3.0/4.0),))
    edge_T3_outter_bottom = T3Instance.edges.findAt(((0,-abs(T3_raduis),tube_len*1.0/4.0),),((0,-abs(T3_raduis),tube_len*3.0/4.0),))
        
    #============移动部件=====================

    currAssembly.translate(instanceList=('T2-1', ), vector=(0.0, -1*abs(1.0*T1_raduis-T1_pair[1]-T2_raduis), 0.0))
    currAssembly.translate(instanceList=('T3-1', ), vector=(0.0, abs(1.0*T1_raduis-T1_pair[1]-T3_raduis), 0.0))
    currAssembly.translate(instanceList=('rigidPlate-1', ), vector=(0.0, 1.0*T1_raduis, tube_len/2.0))
    currAssembly.translate(instanceList=('rigidPlate-2', ), vector=(0.0, -1.0*T1_raduis, tube_len/2.0))
    
    #=============Tie约束施加，区域要求输入元组===============
    t_min = min([T1_pair[1],T2_pair[1],T3_pair[1]])
    tie_len = 0.02 * t_min # 绑定发生的距离是最小厚度的2%
    currModel.Tie(name = 'tie-T1_inner_top-T3_outter_top', 
            main=tuple(edge_T1_inner_top), secondary=tuple(edge_T3_outter_top), 
            positionToleranceMethod=SPECIFIED, positionTolerance=tie_len, adjust=ON, 
            tieRotations=ON, thickness=ON )
    currModel.Tie(name = 'tie-T2_outter_top-T3_outter_bottom', 
            main=tuple(edge_T2_outter_top), secondary=tuple(edge_T3_outter_bottom), 
            positionToleranceMethod=SPECIFIED, positionTolerance=tie_len, adjust=ON, 
            tieRotations=ON, thickness=ON )
    currModel.Tie(name = 'tie-T1_inner_bottom-T2_outter_bottom', 
            main=tuple(edge_T1_inner_bottom), secondary=tuple(edge_T2_outter_bottom), 
            positionToleranceMethod=SPECIFIED, positionTolerance=tie_len, adjust=ON, 
            tieRotations=ON, thickness=ON )
    # abaqus2024修改key word,master和slave更改为mian和secondary
def rotateTT(mdb_name, degree):
    a = mdb.models[mdb_name].rootAssembly
    a.rotate(instanceList=('T1-1', 'T2-1', 'T3-1'), axisPoint=(0.0, 0.0, 0.0), 
        axisDirection=(0.0, 0.0, 1.0), angle=degree)

def globalMeshSize(T1_pair,T2_pair,T3_pair,tube_len):
    '''
    本函数将决定一个全局网格尺寸
    '''
    # 全局网格划分策略的原则：
    # 1. 壁厚方向至少两层，最小壁厚只能两层
    # 2. 不大于四分之一周长的5 %
    # 3. 不大于长度方向的5 %
    # 4. 尺寸尽量设置为整数，如果不能满足 1. 的条件则多保留一个小数位，依次类推

    t_min = min([T1_pair[1],T2_pair[1],T3_pair[1]])
    len_max = tube_len * 0.05          # 不超过管子长度的5%
    cl_max = min([T2_pair[0],T3_pair[0]]) * 3.14 /4.0 * 0.05    # 不超过最小环四分之一周长的5%
    mesh_size = round(min([len_max, cl_max]))
    i = 0

    while t_min/mesh_size >= 2.333 or t_min/mesh_size <= 1.667:
        mesh_size = round(t_min/2,i)
        if t_min/2 == 0:
            mesh_size = 0.4
        i = i + 1
        if i == 3:
            break
    
    while t_min/mesh_size >= 2.333 or t_min/mesh_size <= 1.667:
        mesh_size = mesh_size * 1.1

    return mesh_size

def creatJob(T1_D,T1_t,T2_D,T2_t,T3_t,tube_len,rect_b,rect_h,degree,model_name,job_name):
    # T1_D-2*T1_t=T2_D+T3_D 不输入T3_D,根据几何关系进行计算

    T3_D = T1_D - 2 * T1_t - T2_D 
    import assembly
    TCNTAssembly((T1_D,T1_t),(T2_D,T2_t),(T3_D,T3_t),tube_len,(rect_b,rect_h),model_name)

    currModel = mdb.models[model_name]
    T1_Part = currModel.parts['T1']
    T2_Part = currModel.parts['T2']
    T3_Part = currModel.parts['T3']
    rigidPlate_Part = currModel.parts['rigidPlate']
    currAssembly = currModel.rootAssembly
    T1_Instance = currAssembly.instances['T1-1']
    T2_Instance = currAssembly.instances['T2-1']
    T3_Instance = currAssembly.instances['T3-1']
    topPlate_Instance = currAssembly.instances['rigidPlate-1']
    bottomPlate_Instance = currAssembly.instances['rigidPlate-2']
    
    rotateTT(model_name,degree)

    import material
    #------------铝合金和低碳钢--------------

    mat01_name = "Steel DIN 2393-ST 37_2"
    mat01 = currModel.Material(mat01_name)

    mat01.Density(table = ((7.861E-9, ), ))
    mat01.Elastic(table = ((200e3 , 0.3), ))
    mat01.Plastic(table=((470, 0),
                    (590, 0.078)))

    mat02_name = "alloy 6063"
    mat02 = currModel.Material(mat02_name)
    mat02.Density(table = ((2.78E-9, ), ))
    mat02.Elastic(table = ((68e3 , 0.33), ))
    mat02.Plastic(table=((160.0, 0.0),
                    (205.0, 0.06)))

    assign_mat = mat02_name 
                    
    import section               
                    
    currModel.HomogeneousSolidSection(name='secSolid', 
        material=assign_mat, thickness=None)

    T1_Part.SectionAssignment(region=(T1_Part.cells,), sectionName='secSolid')
    T2_Part.SectionAssignment(region=(T2_Part.cells,), sectionName='secSolid')
    T3_Part.SectionAssignment(region=(T3_Part.cells,), sectionName='secSolid')

    import step

    compute_time = T1_D * 0.71 / 2000.0

    currModel.ExplicitDynamicsStep(name='Step-1', 
        previous='Initial', timePeriod=compute_time, improvedDtMethod=ON)

    currModel.fieldOutputRequests['F-Output-1'].setValues(numIntervals=200)

    currModel.historyOutputRequests['H-Output-1'].setValues(numIntervals=1000)

    regionDef_topPlate = currAssembly.sets['RF-topPlate']
    currModel.HistoryOutputRequest(name='RF-topPlate', 
        createStepName='Step-1', numIntervals=1000, variables=('U2', 'RF2'),
        region=regionDef_topPlate, sectionPoints=DEFAULT, rebar=EXCLUDE)
        
    regionDef_bottomPlate= currAssembly.sets['RF-bottomPlate']
    currModel.HistoryOutputRequest(name='RF-bottomPlate', 
        createStepName='Step-1', numIntervals=1000, variables=('U2', 'RF2'), 
        region=regionDef_bottomPlate, sectionPoints=DEFAULT, rebar=EXCLUDE)

    import interaction

    currModel.ContactProperty('IntProp-1')
    currModel.interactionProperties['IntProp-1'].TangentialBehavior(
        formulation=PENALTY, directionality=ISOTROPIC, slipRateDependency=OFF, 
        pressureDependency=OFF, temperatureDependency=OFF, dependencies=0, table=((
        0.2, ), ), shearStressLimit=None, maximumElasticSlip=FRACTION, 
        fraction=0.005, elasticSlipStiffness=None)
    currModel.interactionProperties['IntProp-1'].NormalBehavior(
        pressureOverclosure=HARD, allowSeparation=ON, 
        constraintEnforcementMethod=DEFAULT)

    currModel.ContactExp(name='Int-1', createStepName='Initial')
    currModel.interactions['Int-1'].includedPairs.setValuesInStep(
        stepName='Initial', useAllstar=ON)
    currModel.interactions['Int-1'].contactPropertyAssignments.appendInStep(
        stepName='Initial', assignments=((GLOBAL, SELF, 'IntProp-1'), ))

    import load

    fixed_region = regionToolset.Region(referencePoints = bottomPlate_Instance.referencePoints.values())
    currModel.EncastreBC(name='fixed-bottomPlate', createStepName='Initial', 
        region=fixed_region, localCsys=None)
        
    cons_region = regionToolset.Region(referencePoints = topPlate_Instance.referencePoints.values())
    currModel.DisplacementBC(name='freedomCons-topPlate', 
        createStepName='Initial', region=cons_region, u1=SET, u2=UNSET, u3=SET, ur1=SET, 
        ur2=SET, ur3=SET, amplitude=UNSET, distributionType=UNIFORM, fieldName='', 
        localCsys=None)
        
    currModel.VelocityBC(name='velocity-topPlate', 
        createStepName='Step-1', region=cons_region, v1=UNSET, v2=-2000.0, v3=UNSET, 
        vr1=UNSET, vr2=UNSET, vr3=UNSET, amplitude=UNSET, localCsys=None, 
        distributionType=UNIFORM, fieldName='')
        
    import mesh

    #------------这里只进行全局网格布种-------------
    # 单元类型：C3D8R，启用 “stiffness” 选项，缩放系数取0.09

    mesh_size = globalMeshSize((T1_D,T1_t),(T2_D,T2_t),(T3_D,T3_t),tube_len)
    
    partInstances =(T1_Instance, T2_Instance, T3_Instance)

    currAssembly.seedPartInstance(regions=partInstances, size=mesh_size, deviationFactor=0.1, 
        minSizeFactor=0.1)

    elemType1 = mesh.ElemType(elemCode=C3D8R, elemLibrary=EXPLICIT, 
        kinematicSplit=AVERAGE_STRAIN, secondOrderAccuracy=OFF, 
        hourglassControl=STIFFNESS, distortionControl=DEFAULT, 
        displacementHourglass=0.09)
        
    meshRegions = (T1_Instance.cells, T2_Instance.cells, T3_Instance.cells )
    currAssembly.setElementType(regions = meshRegions, elemTypes=(elemType1,))

    currAssembly.generateMesh(regions=meshRegions)

    import job
    #-------------此部分提交任务并进行检查--------------
    # 注意开启双精度模式
    mdb.Job(name=job_name, model=currModel, description='', 
        type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, 
        memory=90, memoryUnits=PERCENTAGE, explicitPrecision=DOUBLE_PLUS_PACK, 
        nodalOutputPrecision=SINGLE, echoPrint=OFF, modelPrint=OFF, 
        contactPrint=OFF, historyPrint=OFF, userSubroutine='', scratch='', 
        resultsFormat=ODB, parallelizationMethodExplicit=DOMAIN, numDomains=6, 
        activateLoadBalancing=False, multiprocessingMode=DEFAULT, numCpus=6)

    mdb.jobs[job_name].submit(consistencyChecking=OFF, datacheckJob=True)

    # ----------------------显示设置------------------------
    session.viewports['Viewport: 1'].assemblyDisplay.setValues(renderStyle=SHADED)
    session.viewports['Viewport: 1'].viewportAnnotationOptions.setValues(triad=ON, 
    title=ON, state=ON, annotations=ON, compass=ON)
    session.viewports['Viewport: 1'].setValues(displayedObject = currAssembly)
